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FOREWORD 

This document is one of two documents that constitute the 
final report for Contract NAS8-28089, "Study of Viscous Mixing 
Plume Flow Field." This study was performed by the Lockheed- 
Huntsville Research &: Engineering Center, Inc., for the National 
Aeronautics & Space Administration, George C. Marshall Space 
Flight Center, Alabama. 

The other document is "Theory of Zone Radiometry," LMSC- 
HREC D306101, January 1973. 

The NASA Contracting Officer's Representatives (COR) for 
this contract are Dr. T. F. Greenwood and Mr. D. C- Seymour, 
S&E-AERO-AT. 
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NOMENCLATURE 

Section 2 


A 

scalar 

A 

vector 

A 

tensor 

ANij 

component of the tensor A in the coordinate 
directions WNi, WNj 

APij 

component of the tensor A in the coordinate 
directions Pi, P j 

AXi 

component_of the vector A in the coordinate 
direction U i 

AY i 

component_of the vector A in the coordinate 
direction V i 

AZij 

component of the_tensor A in the coordinate 
directions Wi, Wj 

A0ij 

CSS 

component of the tensor A in the coordinate 
directions W 0i, W0j 

Bij 

matrix function of G 

B 

second-order tensor 

BA 

antisymmetric, second-order tensor 

BS 

symmetric, second-order tensor 

C,CS 

direction cosines 

DA... 

ess — 

components of D in the WA coordinate directions 

D 

tensor (Nth order) 

E, ES 

eigenvalue 
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EBij (EAij) 

EE 

e. 

1 

— i 
e 


FBi (FAi) 
F 

F F 
GC F 
G 


A, G B 
I or G 


-1 


g ij 

H, HS 

T 

M 

MN 

N 

P 

R 

T 

U 

V 

W 


components of E E in the WBi, WBj (W Ai, W A j) 
coordinate directions 

tensor (second-order) 

indicial symbol for WNi 

indicial symbol for W^i 

components of F F in the WBi (WAi) coordinate 
direction 

vector of an antisymmetric tensor 
vector 

cofactors of G 

general-coordinate stretching function, also called 
the metric tensor 

general-coordinate metric tensors 

inverse of G 

indicial symbol for Gij 

indicial symbol for Glij 

eigenvectors 

unit vector 

transformation matrix 
transformation matrix 
transformation matrix 
general coordinate base vector 
position vector 
second-order tensor 
Cartesian-coordinate unit vector 
orthogonal-coordinate unit vector 
general-coordinate unit vector 
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WA.ffB 

WN 

W0 

X 

X I 
Y 
Z 

ZA, ZB 
ZN 
Z (f) 

7 

0 


general-coordinate base vectors 
general, contra variant- coordinate base vector 
general, covariant-coordinate base vector 
Cartesian-coordinate coordinate distance 
indicial symbol for Z i 

orthogonal-coordinate coordinate distance 
general-coordinate coordinate distance 
general-coordinate coordinate distances 
general, contravariant-coordinate coordinate distance 
general, covariant-coordinate coordinate distance 
alternating tensor (third-order) (see Eq. 2.57) 
angle between R and A 

♦ si:*!}:********* 


Lower case letters identify particular components of scalars, 
vectors and tensors 


Section 3 
A, B, C 
B 
Ci 

D 

E 

_F 

G I 
I 


constants 
Nth moment 

any dependent variable on pages 3-12, 3-13; 
constants elsewhere 

rate-of- strain tensor 

decay rate of turbulent energy 

any variable in the blank, and it becomes the 
fluctuating component of that variable 

inverse of the metric G, (see Section 2) 

unit tensor 
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K 

L, Li, Lij 
M 


N 



P 

R 

T 

U 

V 
W 
X 

XG2 

P 

€ 

7 

0 

9 

V 

9 

a 

T 


turbulent kinetic energy per unit mass 
length scales 

any variable in the blank, and it becomes the 
mean value of that variable 

components of the variable named in blanks 
in WN coordinates 

components of the variable named in the blanks 
in the W 0 coordinates 

pressure 

position vector 

temperature 

Cartesian-coordinate base vector 
velocity 

vorticity vector 

Cartesian-coordinate coordinate distance 

lateral distance between two points on the velocity 
profile 

second coefficient of viscosity 
Kronecker delta (see Eq. (3.13)) 
eddy viscosity 

alternating tensor (third-order) (see Eq. (2.57)) 
time 

viscosity 

kinematic viscosity 
density 

pressure tensor 
shear-stress tensor 
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Section 4 
A 


C, C 
L 


1 


i 

N 

P 

Atm 

U 

V 

Y 


6 

€ 


J U 

P 

T 

Subscripts 

<L 

e, oo 

i 

o 



damping factor 
constants 
length scale 

Van Driest mixing length 
correlation parameter 
pressure 

thickness of the mixing region 

velocity along mixing region 

velocity normal to mixing region 

lateral coordinate 

boundary 

eddy viscosity 

molecular viscosity 

density 

shear stress 

centerline values 

values in an undisturbed flow 

inside 

outside 


GENERAL 
denotes functionality 
denotes time average 
denotes vectors 
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denotes tensors 
absolute values 
the determinant 
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Section 1 

INTRODUCTION AND SUMMARY 

Much has been written concerning what is not known about turbulent mix- 
ing. The important thing is that some things are known. Some experiments 
have been performed and some calculation procedures have been developed. In- 
compressible, parallel stream turbulent mixing has been experimentally studied 
and can be accurately mathematically modeled. Afterburning within a rocket 
plume and an injector element's combustion efficiency have not yet yeilded to 
experimental technology or to analytical endeavor. The turbulence literature 
of free shear layers lies between these two extremes — mostly on the side of the 
first example. 

Turbulent kinetic energy models are now a fashionable representation 
of free shear layers. Their use does offer the potential for a significant ad- 
vance in our understanding of turbulence, but this potential has not yet been 
exploited. Furthermore, much of the importance of these models lies in their 
tensor behavior. The mathematical foundations of such behavior is not well 
documented anywhere. Studies to date have been made only on chemically 
simple flows. 

Much maligned eddy viscosity models are still the only ones which can 
be used on combusting flows. Even these models have not been developed to 
their full extent. For example, empirically determined lateral and stream - 
wise variations could be used in an eddy viscosity formulation. (These vari- 
ations would be similar in form to intermittency factors used for boundary 
layers.) However, actual variations have not been satisfactorily determined 
either by experiments or turbulent kinetic methods. 
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The purpose of this report is to give a basic description of the mathe- 
matical tools and models which are presently used to represent turbulent, free 
shear layers. In addition, several recommendations are made for ways in 
which current modeling techniques can be improved. 
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Section 2 

MATHEMATICAL PRELIMINARIES 

Before the laws which govern turbulent flows can be fully applied, one 
must have a proper appreciation of tensor calculus, certain integral theorems 
and averaging procedures, as well as an understanding of ordinary and partial 
differential equations. The purpose of using tensor nomenclature is to ensure 
that complex terms are properly modeled and to simplify coordinate trans- 
formations; it also serves as a type of shorthand, but this is a trivial conse- 
quence. Integral theorems are needed if several control volumes are to be 
related. Averaging procedures are used to convert microscopic fluid prop- 
erties to macroscopic ones. 

This section presents the mathematical tools which are necessary to 
perform fundamental flow analyses. 

2.1 SPACES AND FIELDS 

In fluid mechanics, one is interested in predicting the fluid properties 
density, pressure, temperature, and composition, and the flow properties of 
velocity and head as a function of position and time. Unfortunately, investi- 
gators of other disciplines of interest to the fluid mechanicist, such as num- 
erical analysis and statistical mechanics, have borrowed words to describe 
phenomena of interest to themselves. Therefore terms like vectors, tensors, 
coordinates, spaces, may have several meanings which are only loosely 
connected. 

Historically, geometry and algebra problems were solved before some 
of the more complicated problems in functional analysis; hence, analogies 
were made and terms "generalized." Since this occurrence has now been 
perpetuated for several score years, the current mathematical literature 
has become very lax in separating the real from the abstract. An individual’s 
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unguided initial readings in these subjects are bewildering because a proper 
appreciation of the use of supposedly well-defined technical terms is lacking . 
This section attempts to provide a reading guide for this subject. 


First, a field is defined as a scalar, vector or tensor-dependent vari- 
able which is evaluated over a region of space (call this "position" so "space' 
can be given a different meaning later) and time. Clearly, position and time 
are independent variables. Position is defined by three scalar components 
measured in three coordinate directions. The coordinates need not be recti- 
linear nor orthogonal. Certainly, the distance between two positions is de- 
fined; this distance is said to be measurable by a "metric." 


The ideas of fields have been generalized by analogy into a description 
of a "function space" which may have any number of coordinates and which 
may relate all, some, or none of these coordinates by a metric. 

To make these definitions more meaningful, let us say that three di- 
alects of vector and tensor calculus could be "spoken;" namely, those that 
apply to a field which is completely specified by geometry, to a field which 
is specified by geometry and time, and to a field which is specified by any 
number of coordinates. These fields will be named: Geometric Fields, Un-r 
steady Fields, and Function Spaces, respectively. In Function Space the 
coordinates may or may not be related by a metric, and furthermore, the 
metric may or may not have geometric meaning. 

The language used to describe each of these fields is defined in follow- 
ing paragraphs and then arranged in a form such that it may, subsequently, 
be used. 


2.2 GEOMETRIC FIELDS 

Physical laws must be independent of the units chosen to represent 
physical quantities; hence, units in equations are carefully accounted for. 
Physical laws must also be independent of the coordinates. There are three 
types of geometric coordinate systems: orthogonal rectilinear, orthogonal 
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curvilinear, and general curvilinear. In each of these systems a specific 
set of rules is used to describe the relationship between flow properties and 
geometrical position. These coordinate systems and these rules are de- 
scribed in this section. 

Certain types of quantities which appear in physical laws must be 
classified, because they transform from one coordinate system to another 
by different rules. These quantities are defined and illustrated with ex- 
amples from fluid mechanics. The first of these is a scalar quantity, which 
requires that magnitude only be specified. Examples are: mass, density, 
energy, temperature, volume, pressure, time. The second of these is a 
vector quantity which requires that a magnitude and a direction be specified. 

In fluid mechanics, three scalar values associated with convenient orthogonal 
directions are usually employed to specify a vector. Examples of vectors 
are: velocity, momentum, acceleration and force. 

There are more complicated quantities which require that a larger 

number of scalar quantities be specified. These quantities are called 

tensors, a very general term which may be applied to quantities of any degree 

of complexity. Thus, a tensor of the "zero-order" ("zero-rank") means that 

3° scalars must be specified — this is our "scalar" for a three-dimensional 

space. A tensor of the first-order means that 3 1 scalars must be specified 

2 

as our "vectors." A second-order tensor requires that 3 , nine scalars be 
specified. A third -order tensor requires that 27 scalar components be speci- 
fied. This is the highest order tensor usually necessary in fluid mechanics. 

A second-order tensor requires that a magnitude and two directions be 
specified. In describing the relationship between force, area, and stress, 
the force-vector was described by a magnitude and direction. To describe 
stress, the direction between the area over which the stress acts as well as 
the force must be specified. Examples of second -order tensors are: stress, 
rate-of-strain, momentum -flux, mass moment of inertia. 
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A geometric field is a distribution of scalar , vector , or tensor quantities 
described by functions of space coordinates. The nomenclature for making 
such a description is now to be established. 

In general, upper case letters denote variable names. Sometimes more 
than one letter is used to name a variable. Lower case letters and numbers 
indicate particular components of a variable. Indicial notation was avoided 
because it is difficult to type and to write on a blackboard. The capital letter 
system introduced herein is designed to provide a minimum of symbol changes 
when equations are solved with a FORTRAN program. 

2.3 COORDINATE SYSTEMS 

Coordinate systems are so basic to the understanding of engineering 
problems that their existence and properties are usually presumed to be 
known to the reader. The next few pages attempt to document this "common" 
knowledge. The nomenclature that is needed for subsequent discussion is 
also defined. From the onset, a distinction is made between the coordinate 
systems used in engineering analysis and those which have more application 
to mathematical analysis. The understanding of three engineering and two 
mathematical coordinate systems is essential. 

2.3.1 Engineering Coordinates 

An engineering coordinate system is specified by defining: 

1 . Unit vectors , 

2. Coordinate names, 

3. Stretching functions. 

Unit vectors and coordinates are named with letters and a number; one num- 
ber for each of the three directions. The three common coordinate systems 
are: 
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Cartesian, denoted by U, X as the unit vector and coordinate, 
respectively ; 

Orthogonal curvilinear, denoted by V, Y, 

General curvilinear, denoted by W, Z. 

All of the coordinate systems used herein are taken to have the same 
origin. A point in space is located by a position vector R. In a Cartesian 
system, 


R = U1 * XI + U2 * X2 + U3 * X3 


( 2 . 1 ) 


In any other type coordinate system, the definition of R becomes somewhat 
abstract because the unit vectors vary between the origin and the point; there- 
fore, basic definitions are established for an incremental departure from a 
point, dR. To accomplish this definition, stretching functions which are the 
relationships between increments of coordinates and distances tangent to the 
coordinate directions must be introduced. If a coordinate is an angle, the 
stretching function converts the increment of angle into an increment of length. 
Now 


dR 


U1 *dXl +U2*dX2+U3 *dX3 


VI * (Hl) 1//Z *dYl +Y2 * (H2) 1//2 *dY2+Y3 * (H3) 1 / 2 *dY3 
W1 *(Gll) 1 / 2 *dZl +W2*(G22) 1 / 2 *dZ2+W3*(G33) 1//2 *dZ3 


( 2 . 2 ) 


for each of the coordinate systems, where the G's and the H' s represent the 
stretching functions. The reason for using the numbers in these terms is 
explained subsequently. Note that stretching functions for the Cartesian 
system are unity. 


The values of the unknown stretching functions are determined by first 
calculating (dR • dR) and then formally transforming from the Cartesian sys- 
tem to the curvilinear system of interest. 


2-5 


LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC-HREC D306102 


Recalling that the dot product of two vectors is the product of their 
magnitudes times the cosine of the angle between them, 

(dR • dR) = (dR) 2 

= (dXl) 2 + (dX2) 2 + (dX3) 2 
= HI (dYl) 2 + H2(dY2) 2 +H3(dY3) 2 

= Gil (dZ l) 2 + cos (1, 2)(G11) 1/2 (G22) (dZl)(dZ2) + 
cos ( 1 , 3 ) (Gil) 1 / 2 (G33) 1 / 2 (dZl) (dZ3) + 
cos (2, 1) (G22) 1//2 (G 1 1 ) 1//2 (dZ2) (dZl) + (G22)(dZ2) 2 + 
cos(2, 3) (G22) 1 / 2 (G33) 1 / 2 (dZ2) (dZ3) + 
cos (3, 1) (G33) 1 / 2 (Gll) 1//Z (dZ3) (dZl) + 

cos (3, 2) (G33) 1 / 2 (G22) 1 / 2 (dZ3) (dZ2) + (G33)(dZ3) 2 (2.3) 

where cos(i, j) means the cosine of the angle between Wi and Wj. 

Suppose that the Cartesian coordinates are expressed in terms of new 
coordinates by the equations: 

XI = XI )Z1, Z2, Z3} 

X2 = X2 {z 1 , Z2, Z3} (2.4) 

X3 = X3 |Z1, Z2, Z3^ 

where the braces denote that Xi is a function of Zl, Z2, Z3. The unit vectors 
and stretching functions are not to be considered as independent variables. 

The formal transformation of coordinates is 

dXi = ^ Mmi dZm (2.5) 
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where 


Mmi 


9Xi 

9Zm 


9X1 

9X1 

9X1 

9 Z 1 

9Z2 

9Z3 

9X2 

9X2 

9X2 

9Z1 

9Z2 

9Z3 

9X3 

9X3 

9X3 

9Z1 

9Z2 

9Z3 


(2.6) 


and the summation is performed on the repeated symbol, m. (Note: neither 
unit vectors nor stretching functions have yet been specified in our trans- 
formation scheme. Note also, that the determinant of the matrix M is the 
Jacobian of the transformation.) 


Now the incremental distance can be determined by 


<^> 2 - EEL (Mmi) (dZm) (Mni) (dZn) 


(2.7) 


where the sums are on m, n, i over 1,2,3. It must be remembered that the 
"Mmi" terms in Eqs. (2.5) and (2.7) are matrices, furthermore since the 

dZ terms in Eq. (2.7) do not contain "i", the M's may be summed first. This 
summation of M's is called the "metric," Gmn. It is a matrix; other proper- 
ties of significance that it has are discussed later. 

These operations are summarized in Fig. 1. 

By comparing terms in the two length expressions Eqs. (2.3) and (2.7) 
observe : 

Gki = Gik = (Gii) 1/2 (Gkk) 1/2 cos(i, k) 

therefore the matrix Gmn is symmetric. The significance of these symmetry 
properties is given in the comments section at the end of this section. 
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The Metric — 



(Mmi) (Mni) 


= Gmn = 


9X1 9X1 
_9Z 1 9Z1 


8 

_a 


9X1 9X1 
9Z2 9Z1 


9X1 9X1 
_9Z3 9Z1 


9X2 9X2 9X3 9X3l [9X1 9X1 9X2 9X2 9X3 9X3 ] [ 9X1 9X1 9X2 9X2 9X3 9X3 

+ azl dzl + 8Zl 8Zl] Ifzi 9Z2 + 9Z1 9Z2 + 9Z1 9Z2J [9Z1 9Z3 + 9Z1 9Z3 T 9Z1 9Z3. 

9X2 9X2 9X3 9X3~ 1 [ 9X1 9X1 , 9X2 9X2 9X3 9X3~ ] [ 9X1 9X1 9X2 9X2 9X3 9Xj 

+ 9Z2 9Z1 + 9Z2 9ZlJ [9Z2 9Z2 + 9Z2 9Z2 + 9Z2 9Z2j \dZ2 9Z3 + 9Z2 9Z3 9Z2 9X3 

9X2 9X2 9X3 9X31 [9X1 9X1 , dX2 9X2 9)£3 9X3~ ] [ 9X1 9X1 9X2 9X2 9X3 9X3 

9Z3 9Z1 + 9Z3 9Z1J [9Z3 9Z2 + 9Z3 9Z2 + 9Z3 9Z2j [9Z3 9Z3 T 9Z3 9Z3 9.Z3 9Z3 


ro 

i 

00 


( Gil 
G21 
G31 


G12 

G22 

G32 


G13\ 
G23 j 
G33/ 


Note: Gnm = Gmn; i.e., matrix is symmetric. 
The summation is on i. 


b. Distance — 

(dR) 2 = (Gil) (dZl) (dZl) + (G12) (dZl) (dZ2) + (G13) (dZl) (dZ3) + 

(G21) (dZ2) (dZl) + (G22) (dZ2) (dZ2) + (G23) (dZ2) (dZ3) + 

(G31) (dZ3) (dZl) + (G32) (dZ3) (dZ2) + (G33) (dZ3) (dZ3) 

c. IlGnmll = G11G22G33 +G12G23G31 +G13G21G32 - G13G22G31 - G11G23G32 - G12G21G33 


_ 9CX1.X2.X3) 
" 9 { Z 1 , Z2, Z3) 


the Jacobian 


Fig. 1 - Properties of the Metric 
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To complete the transformation, unit vectors in the curvilinear co- 
ordinate system must be determined. Eq. (2.2) may be written as 


= ||j dZ1 + |f- 2 dZ2 + Us «Z3 


( 2 . 11 ) 


but 


R = U-l XI + U2 X2 + U3 X3 


( 2 . 12 ) 


Therefore, 


W 


wrin l/2 _ Mi U! +M2t; 2 + M1u3 
1 ( G 1 x ) - 9Z1 u 1 + azi + 3Z1 


(2.13) 


W2(G22) 1/2 = |fi U l +H|u2+|i|u3 


(2.14) 


W3 (G33) 


‘/ 2 . t^U.+f^+fffus 


az3 


8Z3 


(2.15) 


The orthogonal curvilinear system is a special case of the general 
curvilinear system which is developed in these pages. Its development would 
show 

Hi to equal Gii 


and 


Gij = 0 for i / j . 


Transformations could have been made directly between any two co- 
ordinate systems, but it is convenient to consider the Cartesian system as 
the reference. (Be cautioned that all literature in this field does not use 
this reference.) 

2.3.2 Mathematical Coordinates 

There are several other coordinate systems vhicn may be defined in 
terms of the metric Gmn. In general they are non -orthogonal , so they are 

named W , Z where additional letters are supplied to name each specific 

coordinate system. These mathematical coordinate systems use only base 
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vectors and coordinate increments; i.e., the stretching functions are not 
defined explicitly but are part of the base vectors. Two additional systems 
are used: those tangential to coordinate lines — denoted by N in the second 
letter of their name and those normal to coordinate lines — denoted by . 


If base vectors are defined by 

WNi = (Gii) 1//2 Wi (2.16) 

they will be parallel to Wi (and also to the coordinate i) but'will not be of 
unit magnitude. Coordinate lengths will be denoted by ZNi; hence, incre- 
mental distances become 

dR = WNI dZNl + WN2 dZN2 + WN3 dZN3 (2.17) 

Notice that dZNi = dZi and that stretching functions are not used. The base 
vector transformations to a Cartesian coordinate system are: 

— — U - I Q ______ 

WN J = aZNj Ui sum 1 = l > 2 > to 6 et j ec l s • (2.18) 


Hence the WN base vectors can be obtained from the U unit vectors by a 
linear -differential-transformation by using the matrix: 



/ 3X1 

3X2 

3X3 


I 3ZN1 

3ZN1 

3ZN1 

_ 3Xm 

Nkm “ 3ZNk " 

3X1 

3ZN2 

3X2 

3ZN2 

3X3 

3ZN2 

1 

\ 3X1 

3X2 

3X3 


\ 3ZN3 

3ZN3 

3ZN3 

Notice that WNi . WNj = Gij. 





2-10 


LOCKHEED -HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC-HREC D306102 


Also 

( d R) 2 = ( MNmi dzNm ) ( MNni dZNn ) sum m, n, i = 1, 2, 3 

and 

^ ^ MNmi MNni = Gmn sum i = 1, 2, 3. 

To further confound the reader, another set of base vectors which are 
perpendicular to the coordinate surfaces, as opposed to being parallel to the 
coordinate directions, is to be defined. Xi Cartesian coordinates and Yi 
orthogonal coordinates are both parallel to the i coordinate direction and 
perpendicular to the j and k coordinate surfaces. Only when non -orthogonal 
coordinates are used is there a difference in these directions. 


The new system is in terms of W0, Z 0. The direction requirements 
for this perpendicular coordinate system are: 


WN1 • W02 = 0 
WN1 • W03 = 0 
WN2 • W03 = 0 
WN2 • W01 = 0 
WN3 • W01 = 0 
WN3 • W02 = 0 


( 2 . 20 ) 


Let 


and 


dR = * dZ01 + . * dZ02 + * dz03 


3R 


3Z01 


. 7 1 v* £-J fS Lu I j 

3Z02 3Z03 


dXi = T”' 7 T dZ0j sum j = 1, 2, 3; to get i eqs. 

&Z0j 


( 2 . 21 ) 


( 2 . 22 ) 


where 


Xi = Xi {Z01, Z02, Z03} for i = 1, 2, 3 
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Therefore a matrix like Eq. (2.6) and a metric like Eq. (2.8) can be constructed. 
If base vectors W^i are defined by 

m/H- 9 R 
W r 1 = 3Zpi “ 



Sum j = 1, 2, 3; to get i eqs. 


(2.23) 


Now how are the functions chosen so that the WCpi's will be normal to coordinate 
surfaces? 


In determining the WN, ZN coordinate system, it was sufficient to 
evaluate the relative partial derivatives of the coordinate lines or the ele- 
ments of the metric. Therefore, let us consider other properties of this 
matrix. Other more direct methods could be attempted (see comment). The 
inverse of G is: G G is given a tensor symbol; proof of tensor properties 
is given in Ref. 2-1, p. 107. It is evaluated thusly. 


1 . Given 

( Gil 
G21 
G31 


2. The determinate of G is |[Gnmi| and it is Eq. (2.10). 


G 1 2 

G13\ 


G22 

G23 1 

(2.24) 

G32 

G33/ 
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of G are 

(Ref. 2- 

1. P- 49): 



GCF 1 1 

= 

(-D 2 

(G22 G33 

- G23 

G32) 

GCF21 

= 

(-1) 3 

(G 12 G33 

- G32 

G13) 

GCF31 

= 

(-D 4 

(G12 G23 

- G22 

G 1 3) 

GCF12 

= 

(-D 3 

(G21 G33 

- G31 

G23) 

GCF22 

= 

(-D 4 

(Gil G33 

- G31 

G13) 

GCF32 

= 

(-D 5 

(G 1 1 G23 

- G21 

G 13) 

GCF 1 3 

= 

(-1) 4 

(G21 G32 

- G31 

G22) 

GCF23 

= 

(-D 5 

(G 1 1 G32 

- G31 

G 12) 

GCF33 

- 

(-1) 6 

(G 1 1 G22 

- G21 

G 1 2) 


(2.25) 


4. Define 


( GCF11 CGF12 GCF13> 
GCF21 GCF22 GCF23 
GCF31 GCF32 GCF33; 


(2.26) 


5. The transpose of 


GCF 


- 


/ GCF 1 1 
GCF12 
\ GCF 1 3 


GCF21 

GCF22 

GCF23 



(2.27) 


6. Lastly, 



/ GCF 1 1 

GCF21 

GCF31 


/ IlGnmll 

IlGnmll 

IlGnmll 

s - 1 = 

GCF12 

GCF22 

GCF32 

IlGnmil 

II Gnmll 

IlGnmJI 


\ GCF13 

GCF23 

GCF33 


\ IlGnmll 

IlGnmll 

IlGnmll 


\ 


) 


(2.28) 


Let G " 1 = GI 


(2.29) 
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Now if GI is chosen as the metric, then the W<$, Z0 coordinate system will be 
correctly aligned with the normals to the coordinate surfaces. This implies 
that the U, X and WN, ZN coordinate systems must both be considered before 
the W0, Z^ can be defined; i.e., elements of GI are calculated from these 
other two systems. Hence; 

E 3Xi 3Xi _ 

az^J az^k “ 

i 

and 

Gljk = W0j • W^k 

On occasion, mixtures of coordinate systems are used. For example, 
a WN1, ZN1; WN2, ZN2; W03, Z03 system could be used. These possibilities 
are considered again on subsequent pages. 

2.4 THE ORIGIN OF TENSORS 

Discounting zero and first-order tensors which are completely described 
by vector analysis, when do tensor quantities arise in engineering analysis? 
Products of vectors without dot or cross simplifications and division by vectors 
are not defined in vector analysis. These quantities arise in the physical 
world and in the mathematical world. They are second (and higher) order 
tensors. 

The product of two vectors is a second order tensor called a dyadic. 

Other quantities have also been given special names, but they are not reviewed 
in this report. 

Now vectors and tensors are defined; and, for consistency, scalars are 
also named at this point. 


(2.30) 


(2.31) 
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Zero-order tensor (scalar): A = A, i.e., one term. (2.32) 

First-order tensor (vector): A = (AZ 1)(W1) + (AZ2)(W2) + AZ3(W3) (2.33) 

Notice there are three scalar coefficients, each of which is associated with 
a coordinate direction. AZ1, AZ2, AZ3 may represent the velocity components, 
for example. 

Second-order tensor: A = (Wl) (AZ 1 1) (W 1) + (W2) (AZ22) (W2) + (W3) (AZ33) (W3) 

+ (Wl) (AZ 12) (W2) + (Wl) (AZ 13) (W3)+ (W2) (AZ21) (Wl) 

+ (W2) (AZ23) (W3)+ (W3) (AZ31) (Wl)+ (W3) (AZ32) (W2) 

(2.34) 

Notice the one magnitude and two directions associated with each of the nine 
elements of this quantity. This may represent the stress relationships (AZij) 
between the force directions (the first Wi) and surface directions (the second 
Wj). Also, note the natural use of indices shown by lower case letters (i, j, 
k, etc.,) to represent general terms. The order of the terms in the triplet 
of products has no special significance; however, it must be remembered 
(or established by convention) which (Wi) is associated with force and which 
with surface (or any other pair of related quantities). Therefore the second- 
order tensor may be stated: 

A = 2T * AZi j) (Wi) (Wj) for i, j = 1, 2, 3 (2.35) 

th 

Hence, an "N" - order tensor is: 

s ■ EE • E (AZij. . . . m) (Wi) (Wj). . . . (Wm) for i, j, . . . m = 1, 2, 3 

(2.36) 

There are N summations, N unit vectors, and N indices on Aij. . . m. Obviously, 
for N> 3 some of the unit vectors must always appear more than one time; i.e., 
A1123 is a component for N =4. 
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Tensors have now been defined. They have not been specifically de- 
fined by transformation laws, as they often are. Such transformation laws 
are discussed — in due time — and they do indeed hold true; but the definitions 
as stated are those in which the application of mathematical modeling of a 
physical situation is most obvious. The necessity of another defining mech- 
anism is indicated if spaces rather than fields are of interest. 

The literature on tensor calculus is not consistent in its use of terms. 
Perhaps a mathematician would rather say that he has generalized certain 
results. Nevertheless, both the A's and the AZij's are called tensors. A is 
independent of the selection of coordinates; i.e., it is invariant with respect 
to coordinates, whereas AZij's are not. 

th 

Returning to the definition of "N" order tensor, and indicating the 
use of either Wi's or WNi's. 

... (AZij. . . m) (Wi) (Wj) (Wm) = 

(ANij . . . m) (WNi) (WNj) (WNm) (2.37) 

for i, j, . . . m = 1, 2, 3 

Note the corresponding tensor components AZij. . .m and ANij, ,m are not 
equal, although either, with the appropriate unit or base vectors, will suffice 
to define the invariant A. Unit vectors are also base vectors in orthogonal co- 
ordinate systems. 

The ANij. . m are called contravariant components of the tensor A. 

If the W0 vector had been used, 

I= EE - E (A0ij. . m) (W0i) (W0j). . (W0m) (2.38) 

The A0ij. . m would be called the covariant components of the tensor A. 
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Now to really complicate matters one might choose a mixture of base 
vectors. For example, consider a second order tensor. 

S = EE (APij) (Pi) (Pj) fori, j = 1, 2, 3 (2.39) 


where 


Pi = WNi and Pj = W0j 

Components in the WN coordinate system are named contravariant 
components and in W0, they are called covariant. APij is said to be co- 
variant with respect to W^j and contravariant with respect to WNi's. Often 
they are distinguished by using subscripts and superscripts. Since these 
components are related by the metric, Gmn, they can be interchanged at will. 

The other major nomenclature system of interest is the indicial one. 

The indicial system strives to preserve the same symbol for the invariant and 
for its components; it distinguishes between the coordinate systems which can 
be defined by the metric Gmn by the location of the indices. Since the indicial 
system uses super- and subscripts for the distinctions between components, 
performing transformations is often referred to as "raising and lowering the 
indicies." 

The tensor components which are used by the engineer are the AZij ones. 
The units of these components are the same in all coordinate systems. The 
other coordinate systems, hence the other components, do not have this prop- 
erty. Historically, in attempting to discuss the components of tensors, with- 
out first specifying the unit (or base) vectors which are being used, special 
names have been given to these components. AZij's are called physical tensors 
(Ref. 2~2, real tensors (Ref. 2-3) and true tensors (Ref. 2-4). In fact, Ref. 2 _ 4 
refers to other tensors as pseudo -tensors. Modern usage gives pseudo-tensors 
a slightly different connotation (Ref. 2*5). The preferred terminology for our 
purposes is to refer to A as a tensor and to the components as being physical, 
covariant, contravariant or mixed. This is the system of Ref. 2*5. 
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2.4.1 Scalar Components of Tensors 

There are two ways that the scalar components of tensors can be used. 
The first way is simply to recognize from the form of tensor involved how 
many scalar components are necessary and when that number is specified, 
stop. If a coordinate change is subsequently desirable, these scalar com- 
ponents can be reevaluated in the new coordinate system. There should be 
a red flag flying at this point for, if scalars are invariant, why do they have 
different values in different coordinate systems? The reason is that they are 
not proper scalar invariants but they are scalar components. This point is 
now illustrated for first- and second-order tensors. This distinction between 
scalar invariants and components can be made formally, and the resulting 
scalar invariants can be used to define the tensor. Usually, this is more 
trouble than it is worth because the scalar components are perfectly adequate 
to model transport phenomena problems. However, the aforementioned il- 
lustration is presented to clarify this point. 


If, for example, the physical components of a vector (first-order tensor) 
are known in a Cartesian and in a cylindrical coordinate system to be 


A = (AX 1) (U 1) + (AX2)(U?) + (AX3)(U3) = 
(AY1) (VI) + (AY2) (V2) + (AY3) (V3) 


(2.40) 


where for this illustration, 


VI is the radial direction 
V2 is the angular direction and, 

V3 is the axial direction. 

Now AX1, AX2, AX3, AY1, AY2, AY3 are scalar components of the vector but 
they certainly are functions of the coordinate systems used. What then are the 
scalar invariants ? 
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2.4.2 The Scalar Invariants of a First Order Tensor 

In any coordinate system, the distance between the origin and any 
point is called the position vector, R. If a vector A is defined at each point 
in space, it is said to be a function of R. 

A = A | R 1 (2.41) 

Since A and Rare vectors, they are invariant with respect to coordinate sys- 
tems. 


Note that A and R lie in a plane, the direction of which is defined by 
the cross product of A and R. The angle between A and R lies in this plane. 

To simplify these relationships, consider a "natural coordinate system" con- 
sisting of a cylindrical coordinate system where the position vector (also 
called the radius vector) is taken as the radial coordinate and also the angular 
coordinate origin line. The axial coordinate is taken to be parallel to the nor- 
mal to the plane containing A and R. 

Two scalar invariants are immediately obvious: (1) the magnitude of 
A, | A |, and (2) the angle between A and R, 0. 

8 = cos" 1 [A • R/|A||R| ] (2.42) 

The third scalar invariant must be a number which fixes the direction 
of the plane containing A and R. There is no obvious choice for such a scalar 
invariant. It is tempting to use the magnitude of the cross product of A and 
R as the third invariant, but it can be calculated from the other two; therefore, 
it does not supply an independent relationship. 

Inability to calculate the third invariant means that the axial direction 
of the "natural coordinate system" cannot be determined from A and R alone. 
Perhaps the direction of the cross product of A and R may be thought of as 
the third invariant, but this still does not yield the necessary third equation, 
unless this direction is assigned an arbitrary value — like zero. Therefore 
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three scalar invariants are not available for a simple algebraic transforma- 
tion of the components A from one coordinate system to another. The formal 
linear transformation for components of A can be accomplished, but there is 
a more simple way. 

If components of A and R are known in one coordinate system (say the 
U, X system) and one wishes to determine the components of A in, for example, 
the W, Z system, the following relationships can be used. 

(A • Wl) - AZ1 

(A * W2) = AZ2 (2.43) 

(A • W3) - AZ3 

Of course one must be able to determine W, Z from U, X, but transformation 

laws have already been presented for this determination. 

2.4.3 The Scalar Invariants of a Second-Order Tensor 


Since the components of a tensor can be defined only after the unit or 
base vectors of a coordinate system have been chosen, a coordinate system 
with three unit vectors: Wl, W2, W3, is chosen here. 


Let B represent the second order tensor, that is being discussed. It is 


; system by its "ij" components. 


B 1 1 

B 12 

B 13 ^ 



B21 

B22 

B23 

1 

(2.44) 

B31 

B32 

B33 J 



repre 

sent the 

scalar invariants 

of this tensor is to 


write it as three vectors and then use the vector invariants. To do this, assume 
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that the dot product of B and an arbitrary unit vector H forms a vector which 
equals a scalar, E, times H. 


B • H = E H 

Let I be the unit tensor: 


I 



0 

1 

0 



in the W coordinate system. 


Therefore, 

(I - E I) • H = 0 


(2.45) 


(2.46) 


(2.47) 


This is a classic eigenvalue problem. The solution of this problem would 
give three eigenvalues: El, E2, E3 and three eigenvectors: 

HI = (Cll) (Wl) + (C 12) (W2) + ( C 1 3 ) ( W 3 ) 

H2 = (C2 1) (W 1) + (C22)(W2) + (C23) (W3) (2.48) 

H3 = ( C 3 1 ) ( W 1 ) + (C32)(W2) + (C33) (W3) 

where Cli, C2i, C3i are the direction cosines of H with respect to W. 

Thus 


El = El HI 

E2 = E2 H2 (2.49) 

E3 = E3 H3 
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are the three vector components of B; however, two of the eigenvalues may 
be complex, making two of the eigenvectors also complex. To avoid this 
problem the use of symmetric and antisymmetric tensors is introduced. Note 
that nine and only nine scalar invariants of B can be calculated — be they real 
or complex. Many other scalar quantities can be calculated from these nine 
scalar invariants; in fact any of these sets of quantities can be designated the 
scalar invariants of the tensor B. 

Any second-order tensor can be represented as the sum of a symmetric 
tensor (BS) and an antisymmetric tensor (BA). 

B = BS + BA (2.50) 

where the components of BS and BA in a W coordinate system are 

BSij = (1/2) (Bij + Bji) 

BAij = (1/2) (Bij - Bji) (2.51) 

Repeating the contraction scheme used on the general tensor B to re- 
duce it to vector components, let 

BS • HS = ES HS 
or 

(BS - ES I) • HS = 0 (2.52) 

The solution of this equation yeilds three eigenvalues (ESI, ES2, ES3) and 
three eigenvectors 

HS 1 = (CS11)(W1) + (CS 12) (W2) + (CS13) (W3) 

HS2 = (CS2 1) (Wl) + (CS22) (W2) + (CS23) (W3) (2.53) 

HS3 = (CS3 1) (Wl) + (CS32) (W2) + (CS33) (W3) 
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Hence 


ESI = ESI HS1 

ES2 = ES2 HS2 (2.54) 

ES3 = ES3 HS3 

(Note: If a vector contains a number in its name, there must be a defined 
unit vector which corresponds to the vector; i.e., the index on ESi corres- 
ponds to the ith unit vector (HSi) not to the ith component of the coordinate 
vector, W.) Since the tensor BS is symmetric all of the eigenvalues ESi's 
and eigenvectors HSi's will be real (Ref. 2.3, p. 61). There are three scalar 
invariants for each of the ESi vectors. 

Consider now the antisymmetric tensor. 

BA = BA11 W1W1 + BA12 W1W2 + BA13 W1W3 + 

BA21 W2W1 + BA22 W2W2 + BA23 W2W3 + (2.55) 

BA3 1 W3W1 + BA32 W3W2 + BA33 W3W3 

in general, but BA is antisymmetric. 

Therefore, by definition 


BA11 = BA22 = BA33 = 0 
BA 1 2 = -BA 2 1 
BA 13 = -BA31 

BA23 = -BA32 (2.56) 
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Since BA has only three components, it may be interpreted as a vector 
(Ref. 2-6, p. 24). 

Let 

BA12 = F3 
BA23 = FI 
BA3 1 = F2 . 

or 

3 

Fk €i i k BAi J (2-57) 

i. j= 1 

where 

£ijk = +1, if i , j , k = 1,2,3 

2.3.1 

3.1.2 

= -1, if i, j,-k = 3, 2, 1 
2, 1, 3 

1.3.2 

— 0, if any i, j, k are equal 

where 

F = FI W1 + F2 W2 + F3 W3 (2.58) 

is called the vector of the antisymmetric tensor. Again F has three scalar 
invariants which were previously discussed. 

Two important operations result: 

BA • A = A x F and A • BA = F x A (2.59) 

where A is an arbitrary vector. 

2-24 


LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC-HREC D306102 


Therefore, the eigenvalue problem for the antisymmetric tensor becomes; 


BA • A = EH = A x F (2.60) 


' Therefore 


E = O, i [f.I, -i I FI, where 
i is the vtf. Ref. (2.3), p. 60. 

The most common usage of the split on B is when the derivative of the 
velocity with respect to the position vector is divided into the symmetric 
rate -of-strain tensor and the antisymmetric rotation tensor. 

2.4.4 Transformation of Tensor Components 

If the components of a tensor are known in one coordinate system, they 
can be transformed into components in a second coordinate system by the follow- 
ing relationships. 

Let one general physical curvilinear coordinate system be designated 
by WA, ZA, and GA and a second by WB, ZB, and GB. Then for a first-order 
tensor : 

F F = FA1 * WA1 + FA2 * WA2 + FA3 * WA3 
or 

F F = FBI * WB1 + FB2 * WB2 + FB3 * WB3 (2.61) 


Where 


FBI 


l/2 

/ GB11 \ ' 9ZB1 
\GAli; 9ZA1 



FA2 



(2.62) 
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and FB2, FB3 are given by like relationships, or in general (Ref. 2.2): 


FBi 



(2.63) 


where the summation is on j and an equation is obtained for each i. 


Second-order tensor: 


E E 


-EE 


EAmn * WAm * WAn 


(2.64) 


where the first sum is on m and the second is on n over the range 1, 2, 3. 


EBij = £E 


GBmm 

GBnn 

1/2 

/ 9 Z Bm 

9 Z Bn \ 

GAii 

GAjj 


^ 9Z Ai 

9 ZA j ) 


EAmn 


(2.65) 


where the first sum is on m and the second is on n over the range 1, 2, 3. 



-order tensor: 


D 


-EE . . . DAmn . . . s 

WAn* . . . * WAs 


* WAm * 


( 2 . 66 ) 


DBij . . . i 


EE E 


GBmm 
GAi i 


GBnn . . . 
GAjj 


GBss 

GAff 


1/2 


/ 9ZBm 9ZBn 

\ 9 Z A i 9Z Aj 


9ZBs \ 
9 ZAf ) 


DAmn . . s 


(2.67) 


If only WN, ZN or only W^, Z $ coordinates are used, the ratios of G's 
do not appear in these transformation equations because they have been ab- 
sorbed in the definitions of the base vectors. If Cartesian coordinates are 
used, the G's do not appear because they are unity. 
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These transformation laws can also be used to go to or from WN, 
ZN — ► W0, Z0. The G's will appear in that case. 

The operator, ^ 


V ~ 3{ )/9R. 


Before the vector and tensor quantities that have just been defined can 
be used effectively, their derivatives must also be evaluated. 


This may be done thusly: 


9 A 
9 Zj 



Bij + 


9 Ai ' 
9 Z j 


Wi 


sum on i, giving j equations. (2.68) 


The significance of this equation is that the scalar coefficients of the vector 
represented by the differentiation are defined. Bij is a matrix which is a 
function of G. 


In a similar manner 


-Ok- = E wi w s + Ti i 

sum on i and j, giving N equations 


.. 9 [Wi Wi] 
9Zn 


9[wi Wj]/9 Zn is a function of G. 


(2.69) 


2.5 CLOSURE 


With the terms defined in these pages, two conclusions are possible. 


1. If the conservation laws can be written in one coordinate 
•’ system — even the Cartesian system — they can be written 
in any other system. All of the unit (or base) vectors can 
be removed from the equations by writing them as a set 
of partial differential equations. 
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2. Before Item 1 can be realized, the proper tensor form of 
every variable of interest must be known. This is a dif- 
ficult requirement if empirical data are used, especially 
for the second order tensors which represent shear-stress 
and Reynolds stress. The Reynolds stress term will be 
considered in detail — and the turbulent kinetic energy 
models will result. 


2.6 COMMENTS 

(This material is not necessary for continuity of this chapter, but is 
necessary for completeness.) 

A. Implications of G being symmetric 

G 12 = G21, G 13 = G31, G23 = G32 


(1) HGnmll = Gil 

G22 

G33 

+ 

G12 

G23 

G31 

+ 

G31 

G 12 

G23 

- 

G31 

G22 

G31 

- 

-Gil 

G23 

G23 

- 

G12 

G12 

G33 


=G1 1 

G22 

G33 

+ 

2 

G 12 

G23 

G31 


-(G31) 2 G22 - Gil (G23) 2 -G33(G12) 2 

(2) GCF11 = (G22 G33 - [G23] 2 ) 


GCF21 = 

-(G12 

G33 

- G23 G31) 

GCF31 = 

(G 12 

G23 

- G22 G31) 

GFC12 = 

-(G 12 

G33 

- G31 G23) = GCF21 

GCF22 = 

(Gil 

G33 

- [G3l] 2 ) 

GCF32 = 

-(G 1 1 

G23 

- G12 G31) 
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GCF13 = (G12 G23 - G31 G22) = GCF31 
GCF23 = -(Gil G23 - G31 G12) = GCF32 
GCF33 = (Gil G22 - [G12] 2 ) 

(3) Aside -determine what GI elements are — 

GCF11 G22 G33 -fG23] 2 ^ 1 

HGnmH ^ G11 (g2 2 G33 - [G23] 2 ) + 2 G12G23G31- Gil 

-[(G31) 2 G22 + (G33) (G12) 2 ] j 

This inequality would become an equality for orthogonal coordinates. 

B. The linear differential transformations were used to relate two different 
coordinate systems in the previous discussion. But could different trans- 
formations be used for the W0, Z0 system? 

For example, usually 


Zl 

= Zl 

jxi, 

X2, 

X3j 

Z2 

= Z2 

Jxi, 

X2, 

X3j 

Z 3 

= Z 3 

jxi, 

X2, 

X3j are known, 


Then 


dZi 


3 Zi 
3Xj 


dXj 


sum on j, giving i equations. 


Now if this procedure were applied to the calculation of explicit values of the 

r) R 

U's, terms involving would have to be evaluated; two problems arise. 

U AlN 
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(1) An integral expression for R in terms of W, Z would have to be determined; 
for example, in cylindrical coordinates, R = VR * YR + VZ * YZ. The lack of an 
integral expression for R in terms of stretching functions and W's and Z's is a 
discouraging fact, (2) Simple expressions like 

Ui = ^Xi sum on J> giving i equations 

would not result, because the Wj's are functions of Xj's; they are not constants. 

C. Evaluation of 3( )/SRfor vectors and second-order tensors. 

(1) Definitions and nomenclature from Ref. (2-1). 
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fr 

z 


11 

f9g il 

♦ 

3 g. •* 

. 3j 


_8X J 

3X 1 

3X 1 . 


‘3g i2 

+ Ihi 

3 g- • 
- -ILL 


.. 3X j 

3X 1 

3X 2 . 


■3g i3 

+ !£ii 

- 


3XJ 

3X 1 

3 X 3 ./ 

4' j ) 

e 2 + f 

(*’ j) F 3 


+ g 


12 


+ g 


13 


I e l +• 


i. j = 1, 2; 2, 3; 3, 1; 
1. 3; 2, 1; 3, 2; 
1, 1; 2, 2; 3, 3; 


(2) Introduction of general unit vectors 


a(wi §ii l/2 ) 


Wi 


3X' 


3g 


■J 


1/2 


f(g ijm } e m = f (gHWnog mm l/2 ) 

1/2 _ 


11 


ax J 


+ g. . 
6 11 


= f(g) g 1/2 Wm 

3X J 


a Wi .dtn(g )l/2 l/ 2 

- Wi + f(g) g 


3X 


•J 


3X 


■J 


mm 


Wm 


where f(g) represents a function of g or g^. 
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^ -E 

3X J i = l 


= A1 


Ai 

awi + 

Wi 

3Ai 



ax^ 


ax^ 


awi 


3W2 

+ A3 

aw3 

ax-* 

* A2 

ax^ 

ax^ 


a ai 

+ W2 

a A2 

+ W3 

a A3 

ax j 


ax^ 


ax j 


But 


= (Bljj wi + (B2jj W2 + ( B3 j) W3 


(see Item (3)). B is a matrix function of g_. 


dA 1 Bij + \ W1 + ( B2j + W2 + 


ax J 


9X 


axJ 


B3j + -^-1 W3 

3X J 


(3) Derivatives of unit vectors 


3Wi_ s ,/ jfe <gji> l/Z A Wi+i 

ax J \ ax J / 2 


n, a/2/3g u 


H 


v3X J 
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OW 1 OWt DW [ 

ax 1 ax 2 ax 3 


awi 

ax J 


aw2 aw2 aw2 

ax 1 ax 2 ax 3 


aw3 aw3 aw3 

ax 1 ax 2 ax 3 


(4) Extension to second-order tensors. 


9 (T) _ 3 (Tij Wi Wj ) 

ax n ax 11 


3Tij 

ax n 


[Wi Wj 


+ 


Tij 


awi wj 
ax n 


9T vj \ / W i Wj ) + Tij 


ax 


n 


Wi 9Wj + Wj 


ax 


n 


ax 


where n is a particular coordinate and i, j are both summed. 
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Section 3 

TURBULENT MIXING VIA TURBULENT KINETIC ENERGY 

If the conservation laws for laminar flows are assumed to be true at 
each instant in a turbulent flow, formal averaging of these equations can be 
performed. These operations produce the Reynolds stress-tensor . This 
section attempts to represent the best possible analysis of these averaged 
equations before resorting to empirical data fits. Such efforts make exten- 
sive use of the turbulent kinetic energy, and have come to be called turbulent 
kinetic energy models to distinguish them from eddy viscosity models. These 
methods are still being researched; hence, a definitive technology does not 
yet exist. Simple cases will be developed and general cases outlined. 

The starting place for analysis is the conservation laws. Immediately 
a nomenclature problem arises. No universally accepted system exists. 

Facility to read indicial systems must be established. Therefore a system 
which is flexible enough to serve as a standard and which preserves the sum- 
mation feature of the indicial system will be established and used throughout. 

An extensive listing of symbols is given in the Nomenclature, page iv, to 
facilitate comparisons to standard references. 

3.1 THE CONSERVATION LAWS 

Consider the following statement of the continuity and momentum equa- 
tions for a single component fluid with no body forces acting on it. A Cartesian 
coordinate system U, X will be used; these coordinates are not moving, i.e., 
they are inertially fixed. 


Continuity equation: 


l£ 

ae 


3 

E 

i=i 


9p » Vi 

axi 


= 0 
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p is density, 9 is time, and V is velocity. Notice, our system of using cap- 
ital letters has been compromised, because Greek letters have also been 
used to name variables. This was done to conserve space in writing the 
equation. FORTRAN names (RH0, THET) could have been used and must be 
in programming, but they are read with the same name. 


Momentum equation: 



E 

j=i 


3 Tji 

axj 


3P 

axi 


(3.2) 


Tji are the components of the shear stress tensor, T. One would like to ex- 
press T as a function of the velocity field. This functionality is denoted by: 


T = T 




(3.3) 


Since T is a second order tensor and Y is a vector; this functionality is not 
straightforward. Consider 


av 

9R 


a second-order tensor 


(3.4) 


This tensor may be written in terms of a symmetric and an antisymmetric 
part. Components of these tensors are: 

a vi i (d vi avj\ i (a vi avj\ e 

axj ■ 2 \axj + axi/ 2 \axj ' 9X ij 1 


It has been previously established that an antisymmetric tensor component 
is related to a vector such that: 


Wk 



9Vj \ 

axi / 


(3.6) 
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W represents the rotation of a rigid element of fluid; there is no shear asso- 
ciated with such a rotation,, 


Let 


Dij 



(3.7) 


This deformation or rate -of-strain tensor is related to shear stress, 


by 


Let 



(3.8) 


T 


A I + B D + C D • D 


(3.9) 


for a Newtonian fluid where 


= -(|^ - p ) 


jLt - 0 V * V, B = + jU, C = 0 


(3.10) 


and 


T = 


-{¥ - " B ) 


v • v 


I + fJ D 


(3.11) 


Assuming 0 =0, 


= [- ^ v • V ] I + fJ.D or 

3 1 

1 avif I 

+ fl Dij 




_ - , 9 Vk ... 

[' 3 £i 8 xk 1J 


(3.12) 
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where 6ij is the Kronecker delta and its properties are: 


6ij = 1 if i = j 
= 0 if i/j 


(3.13) 


Notice Tij = Tji; i.e., T is summetric. Pressure may be considered an 
additional stressing force and a pressure tensor is defined as 


o ij = Tij - P 6ij = -P 6 ij + j <jl [Dij - \ ^ 


9 Yk 
3 C-J axk 
k = l 


6ij 


Finally, after noting that i and j may be interchanged 


“axfe = axm ( M D:ii) “ 3 axm(^ £ 

\ k=l 


avk .... 
axk 6jl 


(3.14) 


(3.15) 


The momentum equation becomes, 


3 \ 3 
, avi , v 9 vi \ a , 

P I ae + 2-j Vj axj ) 2^ axj ^ 


j=i 


J=i 


Dji 


2 avk 

3 axk 
k= 1 


6ji 


ap 

"ax! 


(3.16) 


This is the Navier-Stokes equation (without body forces) as given in Hinze 
(Ref. 3-1, p.16) and Schlichting (Ref. 3-2, p.6l). 


Consider the case of constant viscosity. (Details Hinze (Ref. 3-1, p. 17). 




3 


j=l 


2 3 

9 Vi , II 9 avk 

axj axj 3 ax. axk 

1 k=i 


9 P (3.17) 

axi 
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To average Eqs.(3.1) and (3.17), define mean and fluctuating components 


V = VM + VF, P = PM + PF, p = pM + pF. (3.18) 


Average and combine said equations. 


pM 


( 9VMi 

ae 


VMj 

j=l 


avMfl 
axj . 


/ 3 2 

a pm . / a vMi . 

axi ^ \ axjaxj 

\ j=i 



a avMk \ 

"axi ZLi axk / 

k=l / 


_a 

9Xj 


^<pF VFj > VMi j 


(Ejlj pM<VFiVFj>+ -^< P FVFi> 

V j=i 

+ g|j (<pF VFi> VMj) + gfj (<pF VFi VFj>) 


(3.19) 


where < > denotes the time-average of the included products. Note, values 

for < VFi VFj}, <pF VFi},<pF VFi VFj > must be obtained to solve the momen- 

i 

turn equations. Considering the incompressible continuity equation 



(3.20) 


(3.21) 
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9 2 VMi 

axj'a'Xj + 

(3.22) 

A mean body force, FMi, may be added to either Eq. (3.16) or Eq. 
(3.22) if it is needed. Equation (3.22) with the body force is attributed to 
Reynolds and the term < VFi VFj y is called the Reynolds stress. 

Reynolds stress is a second-order tensor, with three of its terms 
appearing in each momentum equation. The next step in the analysis of 
turbulent flows is to represent this tensor in a tractable form. 

3.2 THE REYNOLDS STRESS TENSOR 


The incompressible momentum equation is: 


3 

9 VMi . 9VM: 

___ + ^ VMj 

i =1 


9 VMi \ _ 9 PM , 

“ “ 9Xi M lu 
i=l 


+ E dj (-P<VFiVFj>) 


Before more generality is introduced, the analysis of incompressible, 
constant viscosity flows will be extended. One would wish to express <VFi 
VFj y as fjVMj. Let us denote what the requirements for such a representa- 
tion would mean. Realize that <VFi VFj> may be considered components 
of a second-order, symmetric tensor. 


1 . 


The Reynolds stress tensor cannot be proportional simply to 
the velocity of the fluid, because the equations of motion must 
be invariant under a Galilean transformation. The same tur- 
bulent shear stress must result with a fixed body and moving 
fluid, or vice versa. Hence, <VFi VFj > may be proportional to 


9 VM _ 9 2 VM 
9Xi or 9Xi 9Xj 


(3.23) 


See Monin and Yaglom (Ref. 3-3, p 371). 
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2. For isotropic turbulent flows, the extra stress term 
must behave as an additional pressure term only. 

(See Ref. 3-3, p 388). 

3. VFi VFj^> is a symmetric second-order tensor; 
hence whatever model is chosen to represent it 
must have the same properties as such' a tensor. 

4. The additional terms which appear in the time -averaged 
continuity and momentum equations are single-point cor- 
relation functions, and as such must obey certain sta- 
tistical requirements. 


These properties are satisfied by the following mathematical models 
(Ref. 3-3, p 388). 

If the fluid is moving as turbulent slug flow; i. e., there are no veloc- 
ity gradients, the turbulence will be isotropic and should behave as an addi- 
tion to the thermodynamic pressure. For the Reynolds stress tensor to be 
isotropic, its diagonal elements must be constants and its off-diagonal 
elements zero. To represent this behavior; 

1 3 ? 

p < VFi VFj ) = j <VFkVFk> 6ij = | p K 6ij (3.24) 

k=l 


where K is the turbulent kinetic energy per unit mass (abbreviated TKE). If 
the flow has velocity gradients, < VFi VFj> cannot be simply set equal to a 
constant times the velocity because each term in an equation must have the 
same tensor character. But 

/ 9VMi 9 VM.j 

l 9Xj 9Xi 



which equals Dij, is a suitable second-order tensor; therefore, let 


< VFi VFj > f 


F jvM 


/ 9VMi 9 VMi 

l 9Xj + 9Xi 


(3.25) 
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where F |VMj is a function of the vector VM. Since the isotropic limit pre- 
viously given can be added to Eq. •(3.24) without destroying the functional 
character of the equation 

<VFi VFj> = | K 6ij + FjVMj(^^ + (3.26) 

Fj VM j may be a scalar, a second-order tensor, or a fourth-order tensor. 

It cannot be an odd-order tensor because proper contractions cannot be made 
to maintain the tensor nature of the equations. Even-order tensors of order 
higher than four when contracted would become identical with second- and 
fourth-order terms. 

What has been discussed are possible values of <VFi VFj>. Presently 
an analytical method for obtaining these values will be discussed. Empirical 
determinations of turbulent mixing have been made which fit some of these 
functional forms. These may be stated as follows. 

Fisa scala r. 

<VFi VFj > = ^ 6ij - L K 1//Z Dij (3.27) 

where L is the mixing length. Usually the first term on the RHS may be as- 
sumed negligible (Ref. 3-4, p. 33), then the expression becomes of the form 
of that suggested by Boussinesq. 

F is a second-order tensor. 

3 

p < VFi VFj > = | Kp 6ij - jr (p KV /2 Dkj + Ljk Dki) (3.28) 

' ' k=l 


This form was suggested by Monin. Lij are components of a symmetric 
second-order tensor; hence, they represent six length "scales," since their 
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dimensions are in length. The introduction of Lij is useful to interpret the 
physical phenomena of mixing, but it does not simplify the problem of deter- 
mining the Reynolds stress tensor, either the six components of L or K. YF 
VF> must be known. 

F is a fourth-order tensor. 

3 

p <VFi VFj > = y K p 6ij - p ^ KKijkm Dkm (3.29) 

k, m = l 

where KK must be of the form: 

KKijkm - j ^Kik 6jm + Kjk 6imj (3.30) 

No one appears to have used this as an empirical form of eddy transport. 

In all of these forms, the approximation of neglecting the first term 
on RHS of the Reynolds stress tensor is frequently encountered. 

Now we shall return to the problem of obtaining the Reynolds stress 
tensor by analysis. 

3.3 TURBULENT KINETIC ENERGY 

To obtain expressions for <VFi VFjX the conservation laws already 
described may be manipulated until partial differential equations for the com- 
ponents of interest are obtained; these are called dynamic equations . A more 
satisfying development which results in the same equations may be made by 
using the statistical properties of the variables of interest; specifically, equa- 
tions for the moments are determined. Both of these developments will be 
reviewed. 
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3.3.1 The Dynamic Equations (Hinze, Ref. 3-1, p 250) 


The incompressible momentum equation written in terms of the mean 
and fluctuating terms is: 


3 VMi 
36 


3 VFi 
30 


+ 


E 


k=l 


1 3PF 
p 3Xi 


(VMk + VFk) 


3 (VMi + VFi) _ 2. 9 PM 

3Xk " - p 3Xi 



(vMi + VFi 
3Xk3Xk 


(3.31) 


Assume steady mean motion, = 0, and add the fluctuating part of the 

d o 

continuity equation, Eq. (3.21), to the momentum equation. 



3 VMi 
3Xk 


+ VFk 


3 VMi 
3Xk 


+ VMk 


3 VFi 
3Xk 


+ 


VFk 


3 VFi 
3Xk 


+ 


VFi 


3 VFk 
3Xk 


RHS of Eq. (3.22) 


(3.32) 


Combine the last two terms on the LHS and substract the time -average of 
this equation, Eq. (3.32), to obtain: 


3 VFi 
39 


3 


+E 


^VFk 


3 VMi 
3Xk 


+ 


- <VFi VFk> 


1 3PF 

p 3Xi 


VMk 


3 VFi 
3Xk 


+ — 
3Xk 


VFi VFk 


+ 


v 



3 2 VFi 
3Xk 3Xk 


(3.33) 


There are three equations corresponding to the three momentum equations. 
Denote another one by j and write it out. 


9 VFj 
30 



k = l 


^VFk 


3 VMi 
3Xk 


+ VMk 


9 VFj 
3Xk 


+ — — 

3Xk 


-<VFj 



I 9PP 

p 3Xj 


+ V 


3 


E 


3 2 VFk 
3Xk 3Xk 


VFj VFk - 


(3.34) 
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Multiply Eq. (3.33) by VFj and Eq. (3.34) by VFi, then add, and average. Re- 
place the naturally occurring average terms on the RHS of the resulting 
equation with their equivalents as listed in Hinze, (Ref. 3-1, p 251) to get: 


+ VMk 


k = l 

a 

axk 


a YMj 


A <VFi VFj > +£ I < VFj VFk> 9 VMi 

< VFi VFj > ] =-jr A- <VFi VFj VFk > 

J 1 


+ < VFi VFk> w 
3^ 

r 

kTl 


?(§|-i< PFVF j> + iki < PFVFi >) + 


+ K pF (w j ♦ 


\ 


. „ v* 9 2 <VFi VFi > 

v Lu axk axk 

k=l 


a VFi 
9Xj 


z , 3 VFi a VFj 
axk axk 


(3.35) 


<VFi VFj y represents the mean value of two velocity components or the ele- 
ments of the symmetric Reynolds stress tensor. There are six possible 
unique combinations of i, j; namely, 1, 1; 2, 2; 3, 3; 1, 2; 2, 3; 3, 1. This 
means that Eq. (3.35) represents six equations. 


A contraction (in the terminology of Cartesian tensor analysis) of 
these six equations is represented by 


_a_ 

30 


E 


i = l 


3 

< VFi VFi > + ^ 
i, k = l 


2 VFi VFk 


aVMi 

axk 


+ 


VMk 


S5k < VFi VF1 > 


3 

-*E 

i, k = l 



VFi VFi 



+ 
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3 

+ *I 

i,k = l 


9 < VFi VFi > 

axk axk 


a VFi a VFi 


axk axk. 


(3.36) 


This equation is obtained by adding the three Eqs. (3.35) for which i, j = 1, 1; 2, 

2; 3, 3. (Note that < VFi VFi> is K, which is two times the turbulent 

i = l 

kinetic energy. ) 


To use either the six Eqs. (3.35) or Eq. (3.36), now that they have 
been derived and all of the terms have been identified, the averaged quanti- 
ties; i.e., those in broken brackets, ^ X may be named as new dependent 
variables. This means that the number of dependent variables may be made 
somewhat less than the number of previously named variables. For example, 
p < VFi VFj > may be called pK, two variables, and that p, VFi and VFj (three 
variables) are not necessary. K is the turbulent kinetic energy. 


Attempts to solve for all of the dependent variables at this point, by 
making empirical correlations when necessary, is called "using a turbulent 
kinetic energy approach" to represent turbulence. Various attempts to accom- 
plish this task will be reviewed, but first an alternate derivation of Eqs. (3.35) 
will be outlined. 

3.3.2 The Moment Equations (Monin and Yaglom, (Ref. 3-3, p 374) 

At best the operations involved in obtaining Eqs. (3.35) are cumbersome, 
they certainly are not obvious. A more direct approach is to calculate the 
mean value of an VFi VFj term directly. This can be done by formulating the 
equations for moments. 


Consider n different fluid dynamic variables (dependent variables). 
Generally, these may be functions of position and time. Denote these variables 
by the symbol Ci. In general, different Ci's may be functions of different points 
m space and time. The " n " th moment of these variables is defined by 
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Rl, 9 j * C2 j R2, 9 | * C3 j R3, 0 


Cn 


jRn, 0^ 


(3.37) 


the Rj points may be coincident or different. 

Reversing the order of averaging and differentiation of Eq. (3.37) with 
respect to time. 


9B ✓ 9C 1 ^ 

"90" = < “30- ' C2 


Cn > + < Cl 


9 C2 ... 
99 


Cn y + 


... + <Cj * ... * (3.38) 

th. 

B has been written as a scaler; it may be an "n" order tensor, if the Ci 
variables are vectors. However, the defining conservation equations are 
already written as partial differential equations, which means that one may 
simply consider the vector components as scalers (if no coordinate change is 
used; i.e., if no vector or tensor component has to be calculated in another 
coordinate system). 


To proceed, one uses Eq. (3.38) to define a "B," the partials with re- 
spect to time that appear on the RHS are obtained directly from conserva- 
tion laws. This is the logic for why the operations in the previous section 
are reasonable. 


The previous section may be summarized by saying that the time - 
averaged values of the single-point correlation functions, <( VFi VFj /, were 
formulated as a set of partial differential equations. The correlations were 
single -point because all of the position vectors in Eq. (3.37) when applied to 
this problem were identical. The determination of an average value may be 
described as determining the correlation of the variables involved. More 
formal definitions of this term can be made (Ref. 3-3, p 228), but these do 
not contribute to our further understanding of turbulent mixing. 
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3.3.3 Closure 

Returning to the problem of obtaining a solution to either Eqs. (3.35) 
or (3.36), one sees that an attempt to calculate directly a correlation function 
leads to difficulties because more correlations are introduced into the equa- 
tions. Thus more variables than equations result — this causes what is termed 
the closure problem . Namely, new equations for the unwanted correlation 
functions must be specified before a solution can be obtained. Reference 3-3, 
p 377 summarizes the history of recognizing such difficulties. 

The next section reviews the better attempts to supply such new equa- 
tions. 

3.4 SOLUTIONS TO THE TKE EQUATIONS 

The dynamic equations Eqs. (3.35) whichrepresent the velocity correla - 
tions, < VFi VFj >, are one starting place for using empirical turbulence data, 
but it is not an obvious starting point. However, Prandtl reportedly suggested 
an equation for turbulent kinetic energy which was "derived by logical reason- 
ing and dimensional analysis," Eckert and Drake (Ref. 3-5, p 369). This equa- 
tion is similar to the contracted form of the general kinetic energy equation, 

Eq. (3.36). In any event, nothing more can be accomplished until some empir- 
ical information is introduced. There are four ways that this can be accom- 
plished: (1) postulate enough information to solve Eqs. (3.35); (2) solve Eq.(3.36) 
or one similar to it for K, using empirical data as necessary; (3) repeat method 
(2) except solve for both K and an additional length scale, L, from equations 
like Eq.(3.36); (4) solve Eqs. (3.35) for the components of the Reynolds stress 
tensor and one equation like Eq. (3.36) for a length scale. Without posing the 
question why, let us examine what these four methods would entail. 


3-14 


LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC-HREC D306102 


3.4.1 Solutions for all Correlations 


In order to understand turbulent kinetic energy methods, consider the 
work reported in Ref. 3-6, which is a somewhat more general account of that 
reported in Refs. 3-7 and 3-8. 


Rearranging Eqs. (3.35) and using Eq. (3.21) 

-^<VFiVFj> + -gffk |(vMk< VFi VFj>j + <VFi VFj VFk>j = 

k ~ X 


3 

£ 

k = l ^ 


< VF i VFk> 2™- 1 + <VFi VFk> '^Mi 


(3.39) 


+ [The last six terms on the RHS of Eqs. (3.35)] 


The first pair of these last six terms can be moved to the LHS of our 
new equation. The second pair are to be left on the RHS. The last pair con- 
tain a term that with the definitions on p 251 of Ref. 3-1 can be written: 



But 


, 9 (VFj Tki) 
9Xk 


VFj 


9 7 ki 

9Xk 


Tki 


9 VFj 
9Xk 


3-15 


LOCKHEED -HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC-HREC D306102 


Therefore all of these operations can be combined to give: 


9< VFi VFj > 
99 


+ 


3 


z 


9 

9Xk 


VMk < YFi VFj > + < VFi VFj VFk> + 


+ VFi (6kj PF - rkj) + VFj (6ki PF - rki) 


3 



< VFj VFk> 


9 VMi 
9Xk 


+ < VFi VFk > 


9 VMj 
9Xk 


+ 




9 VFj 
9Xk 


S)* + <f k j 


9 VFi\ 
9Xk/ 


(3.40) 


Since Eqs. (3.40) have introduced new correlation terms, models must be used 
to relate all of the terms on the RHS to <VFi VFj> and to mean flow properties. 
Mellor and Herring (Ref. 3-6) assumed, on the basis of physical arguments, that 


<VFi rkj)+ < VFj rki > 


( 9 < VFj VFk> 9< VFk VFi> 9<VFi VFj > 
V \ 9Xi 9Xj 9Xk 

(3.41) 



9 VFi 
9Xk 



9 VFj N 
9Xk / 


2 K 3 ^ 

3 L4 


6ij 


(3.42) 


Where L4 is a function of v only near a wall. 


I 

P 



9VFi \ 
9Xj / 



K 1/2 / 

r 7l (<VFiVFj> 



(3.43) 
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<p FV Fi> _ K l/2 L2 8K 
vm; - 2 axi 


(3.44) 


< VFi VFj VFk> = K 1 ^ 2 L3 ( 9 + 


3<VFi VFk> 3<VFi VFi> 

axi axk 


(3.45) 


With these five assumed functions (Eqs.(3.41) through (3.45)), Eq. (3.40) becomes 


D<VFi VFi > 
D0 




H K 1 / 2 L3 V aiV|lVFk> + 


a<VFkVFi> 3<yFi VFi>]| , K 1 ^ 2 L2 / 

axj axk )\ 2 l 

( 5k i Hr + ^ - 

-<VFkVFj> ^iV I- K^<VFiVFj>- 



2 K 3 / 2 

3 L 4 


6ij 


(3.46) 


The LI, L2, L3, L4 are empirical functions which must be specified. 

Notice that Eqs. (3.46) can be contracted; i.e., the three equations for 
i, j = 1,1; 2, 2; 3,3 are added, and divided by two to give; 

^ -t (ak[(^« l/2 -)(^ ♦ 
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+ E 8 <W— )+ k1/2l2 ^# ) ]- 

i = l ' 

_£<VFkYFi> 2J^V ^ (3-47) 

There is not an immediate application for this equation, since it represents 
less information than Eqs. (3.46). Equation (3.47) is known as the turbulent 
kinetic energy equation. 

The six equations represented by Eqs. (3.46) together with the mean 
continuity equation, Eq. (3.21), and the three mean momentum equations, Eqs. 
(3.22), can be solved for the ten unknowns: VM1, VM2, VM3, PM,( VF1 VF1>, 
< VF2 YF2 >, <VF3 VF3>, <VF1 VF2>, <VF2 VF3>, <VF3 VF1>. if the neces- 
sary boundary conditions are specified. This constitutes a solution to the 
closure problem. 

3.4.2 Boundary-Layer Applications 

Generally, the Li's and the necessary boundary conditions are not 
known; hence, further progress is made by restricting the class of flows one 
tries to describe. Remember the restrictions to this point are: (1) steady 

mean flow, (2) constant viscosity, and (3) incompressible. 

The largest subclass of flows which have been investigated is two- 
dimensional thin-shear layers. Notice these flows are not necessarily bound- 
ary layers. For this class of flows, let: VM3 = 0 and the primary flow 
direction be XI. Hence, 


9VM1 9YM2 
9X1 0X2 " 


(3.48) 


3-18 


LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 



L.MSC-HREC D306102 


DVM 1 


1 3PM + 3 


D Q ~ p 3X1 3X2 


^ -<VF1 YF2 > + V j (3.49) 


D<VF i > 3_ 

D0 ” 3X2 




VM1 


K 

3 


\ i Arp 1 vpov 3 VM1 2 K 

) " ^ WF1 VF2 > ax2 “m 


3/2 


(3.50) 


D <VF2^ > 
D0 


3X2 


3 ^ + K 1//2 L 3 ^ 


3<VF2 / ~l/2 ? 3K_ 

3X2 + ^ ^ 3X2 


K 


1/2 


3 LI 


( 


VM2 


2 K \ 2 K 3 / 2 

' 3 / 3 L4 


(3.51) 


D <VF3 2 > 


D0 


3 

3X2 


[( „ ♦ «) ^<f§^] - (VM3^ 


K \ 2 K 


) 


3/2 


3 L 4 


(3.52) 


D<VF 1 VF2> 

D0 3X2 


■ 1/2 


8 >2( V+ K 1/2 L3^ <VFIVF2> 


3X2 


<VF2 2 > dVM^l _ K_/_ <VF1 VF2> 


(3.53) 


Donaldson and Rosenbaum (Ref. 3-9) obtained equations quite similar to these. 

The variations (Ref. 3-9) introduced were: LI = L3, a different L4 and L2, re- 
placing K 3 ^ 2 with < VFi 2 > in Eqs. (3 . 50 and 3 . 52) and added an Y^ 2 ^term 

to Eq. (3.53). These minor differences were due to using slightly different as- 
sumptions to remove the unwanted correlation terms from the dynamic equations. 
Equations 3.48 through 3.53 were solved by Ref. 3-9 for specified initial condi- 
tions and zero correlation coefficients at the wall and at the edge of the boundary 
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layer. Pressure gradients along the boundary layer were also specified. 
Note that the specification of the Li functions and the boundary conditions 
made this a boundary layer analysis instead of a free shear layer analysis, 
i.e., the basic equations were not changed. 


The interesting feature of this boundary layer solution was that prop- 
erties not easily predicted were calculated (Ref. 3-10). Specifically, transi- 
tion was calculated for arbitrary perturbations imposed on an initially laminar 
boundary layer. 


3.4.3 Solutions for the Contracted TKE Equation 


Recognizing that solving Eqs. (3.48) through (3.53) was not a trivial task, 
Mellor and Herring (Ref. 3-9) argued that a nearly isotropic turbulence could 
be defined by a limiting form of Eqs. (3.46) and (3.47), namely: 


< VFi VFj> = 6ij -j 


K 


1/2 


LI 


/aVMj 9 VMi 
l 9Xi 9Xj 


(3.54) 


and 


D(K/2) 

D0 


3 

E 

k = l 


9Xk 


1 ( 5 - 


+ | K 1//2 L3 + K 1//2 L2 


\ 9(K/2) ] 

/ 9Xk J' 


3/2 


N 9 VMi \ K 

<VFkYFi> ^ - -£4 


(3.55) 


This is gratifying, since an equation like Eq.(3.54) has already been postulated 
as a mathematically possible Reynolds stress model. Equations (3.54) and 
(3.55) also supply a satisfactory solution to the closure problem, again provid- 
ing the Li's are specified. 
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The two-dimensional thin shear layer version of Eqs. (3.54) and (3.55) 


is 


<VF 1 VF2> = 


- * 1/2 « *%£ 


(3.56) 


and 


D(K/2) 

D0 


9 

9X2 


+ K 1//2 LI 


^K 1/2 (L2 + + 

(dVMl\ 2 K 3//2 
l 9X2 ) " L4 


(3.57) 


From numerical experiments in Ref. 3-6 


Ll = L2 + j L3 


(3.58) 


Assume V negligible. 


D(K/2) _ _9_ /f K l/2 n 9(K/2 )\ 
D9 ‘ 9X2 \L X J 9X2 / + 

+ k ‘/ 2 L1 - # 


Equations (3.56) and (3.59) are precisely the same as those obtained 
by investigators who postulated physical models for the conservation law con- 
trolling the supply of turbulent kinetic energy (Ref. 3-5, p.366). 

Now expressions for the Li's must be determined. References 3-6 and 
3-9 merely assumed them. Assuming Ll proportional to L4 and recognizing 
that (K 3//2 /L4) represents the decay rate of turbulent energy, E, Rotta (Ref. 
3-8) suggested using an additional transport equation to calculate this 
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decay rate. Rather than use strictly empirical information to describe L4, 
Rotta (Ref. 3-11) defined L4 in terms of a two -point correlation function. Upon 
approximating this correlation function for boundary -layer type flow, a decay 
rate was defined and a transport equation for such a decay rate was written. 

A similar equation was used by Spalding and associates (Ref. 3-12). 


DE 

D0 



9 

axj 


§<VFj VFk> 


9 E 
9Xk 


3 



j»k = l 


(<VFj VFk> | - G3 E 2 /K 


(3.60) 


Where Ci's are empirically determined constants. 

Boundary conditions must be specified for Eq. (3.60) as well as all other 
equations used to represent turbulent transport. Seldom are such conditions 
known, but both Refs. 3-9 and 3-11 indicate that reasonably obvious and simple 
boundary condition estimation leads to valid flowfield representations. Of 
course, better answers can be obtained if exact boundary conditions are known. 

In summary, the analysis of turbulent kinetic energy which was made 
by Rotta (Refs. 3-7 and 3-8) formed the basis of a very powerful tool for de- 
scribing turbulent transport. Unfortunately, this work was not utilized by 
English-reading investigators until Bradshaw's (Ref. 3-13) and Glushko's (Ref. 
3-14) work became known. A convenient review and slight extension of these 
works is Ref. 3-15. Since then three types of turbulence models have evolved: 

1. Equation (3.55) is solved for K for algebraicly specified 
length scales, the Li's. 

2. Equation (3.55) is again solved for K and hence, Eq. (3.60) 
is used to calculate E; no Li's are necessary. 

3. Equations like (3.46) are solved for the components of 
<VFi VFj> with either all Li's specified or Eq. (3.60) being 
used for L4 and the others specified. 
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Models 1 and 2 offer no material advantage over eddy viscosity models, 
although they are not much more computationally difficult and they may be 
somewhat more accurate. 

Model 3 is extremely useful because it offers flexibility to estimate 
certain effects with only a minimum of specified boundary conditions. How- 
ever, general computation schemes utilizing this technique have not been 
extensviely investigated. 

Turbulent kinetic energy models may be used to represent both boundary - 
layers (Ref. 3-16) and free shear layers (Ref. 3-17), provided that the appropriate 
empirical information for either such flow is utilized. A "best" model has not 
yet been established, nor is it likely to be. The type of choices of modeling 
terms is illustrated by the discussion of Harsha (Ref. 18, p. 357). 

It should be remembered the variable density flows and flows for which 
non-Cartesian coordinate systems are desirable have not yet been discussed; 
they will be in subsequent sections. Note also that for the description of in- 
compressible flows, the introduction of the conservation of energy equation is 
not necessary. 

3.5 FLOWS WITH VARIABLE DENSITY AND COMPOSITION 

The procedures which were used to analyze single -component, incom- 
pressible flows can conceptually be extended to handle multicomponent, com- 
pressible flows. The effect of compressibility in a chemically well-behaved 
flow probably will be predictable in the near future; several successful re- 
searches have already been reported in this area. However, the inclusion 
of complex chemical effects in turbulent kinetic energy models is a long way 
from being accomplished, and it is probably not worth the effort to try to do 
so. 
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Density can be described with mean pM and fluctuating, pF, components 
for a single-component fluid. An analysis using either the philosophy of calcula- 
ting the dynamic equations or simply a direct determination of all of the moments 
of interest would result in a set of ten equations; six like those in Eqs. (3.35) 
except for some additional terms on the RHS and four new ones for the var- 
iables <pF VFi)> and <pF pF>* If the thermal equation of state is used, 
temperature, both mean, TM, and fluctuating, TF, may be used to replace 
the density variables. Donaldson et al. , (Ref. 3-19) report a set of dynamic 
equations and then restrict the set to a simple enough atmospheric shear flow 
that a closed set of turbulent correlation equations is obtained. Initial condi- 
tions are postulated and the correlation and conservation equations are solved. 

Laster (Ref. 3-20) derives the dynamic momentum and continuity equa- 
tions, uses a Crocco type energy equation, and restricts the dynamic equations 
to a boundary layer type flow to obtain closure. The empirical models assumed 
cause Laster's model to be hyperbolic — a situation to be avoided if possible. 
Heyman (Ref. 3-21) proposed a somewhat similar model, but one which was 
much easier to use. Heyman admits extensive use of empirical data, but 
nevertheless reports some very impressive data comparisons. 

Another alternative is that the effects of density fluctuations can be 
ignored and then models based on incompressible flows can be used directly. 
References 3-12, 3-22 and 3-23 report good data comparisons using this phi- 
losophy; this is not too surprising. Reference 3-1 p 19, suggested such slight 
dependence based on order of magnitude estimates. 

Reference 3-24 is a most interesting paper which formulates a very 
complete set of dynamic equations. Not only are energy and composition 
fluctuations considered, but also magnetic and electrical. This formulation 
assumes a boundary -layer type flow; hence, the final set of closed equations 
looks like the kinetic energy equation, Eq. (3.36), for which empericial length 
scales are to be supplied. However, they obtained extra equations for thermal 
and mass diffusion; these equations contained parameters analogous to the 
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length scales. They concluded that there were insufficient available data to 
estimate these necessary empirical parameters. 


3.6 FLOW DESCRIPTIONS WITH NON-CARTESIAN COORDINATES 


The real benefit derived from using turbulent kinetic energy models 
is that by properly modeling the Reynolds stress tensor, empirical data taken 
in one flow geometry can be used directly to describe a different flow geometry. 
As the discussion in Section 2 implies this is not a particularly easy task to 
accomplish. Furthermore, very few problems are worked through to their 
solution in general coordinates; i.e., if there are errors in the general co- 
ordinate formulation, they might be "washed out" in a simple problem. Still, 
an attempt at a general coordinate formulation has been reported (Ref. 3-9). 

An important subclass of coordinates is a cylindrical coordinate system 
with axial symmetry. The dynamic equations are appreciably simpler to for- 
mulate in such a system, but, if one’s attention is limited to two-dimensional, 
thin shear layer flows, the physical development of the turbulent kinetic en- 
ergy equation can be easily accomplished (see Refs. 3-12 and 3-18). In fact, 
the general dynamic equations referenced in the previous paragraph have also 
been reduced to those for an axial symmetric coordinate system (Ref. 3-25). 

The turbulent kinetic energy models suggested in Refs. 3-12 and 3-22 are 
valid for plane and axisymmetric flow. 


Some details of the general coordinate formulation are now reviewed 
to indicate the type of difficulties encountered in such studies. The Ref. 3-9 
formulation of the dynamic equations is in indicial rotation. These equations 
are written with the nomenclature of Section 2. They are, for incompressible 
flow: 


3 

p <YF0i VFjZ$k > + p £ VMNj <VF0i VF0k> = 

j=l 
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-p ^ <VFNj VF0k> 
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■p X <VFNj VF ^ i> 
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9V ^ ~ _P X gtj <VFNj VF0i VF0k> 
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avFOki 

3Zi 


+ fj, GImn 

m,n = l 


9 Zm 9 Zn 


<VF0i VF0k> - 2 ^ 


m,n-l 


(3.61) 


Notice : (1) Coordinates Z, ZN, Z0 are taken as being identical; (2) both VFNi 
and VF0i components of vectors are used in these equations; (3) due to sym- 
metery, six equations are represented by Eqs. (3.61); and (4) Eqs. (3.61) are 
in terms of the covariant components of the Reynolds stress tensor. To use 
Eqs. (3.61), one would first rewrite them in terms of the physical components 
of the Reynolds stress tensor. However, since Ref. 3-9 also modeled the cor 
relation terms in the RHS of Eq. (3.61), let us simply study their result — in 
covariant form. Also the continuity and momentum equations are written for 
completeness. 
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j=l 


9 VMj _ 
9Zj - 


(3.62) 
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39 


£ dh ( TN ^ ji " p <VFNj VF ^ i >) 
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<VF0i VF0k> + 2^ VMNj "Sf/ <VF0i VF0k> = 


(3.63) 
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(3.64) 


The significant feature is that all geometric effects have been removed once 
the length scales are determined in any one experiment. Note that the need for 
the introduction of physical components still exists in Eqs. (3.64). 
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3.7 SUMMARY 

The axial symmetric model for thin, two-dimensional shear layers 
is probably the only one that will be used in the near future. Truly complex 
three-dimensional turbulent flows can be modeled with Eqs. (3.62) through 
(3.64), if sufficient theoretical or experimental work is ever accumulated to 
evaluate L4 and L5. A synopsis of this section is given in Table 3-1 by 
listing a set of typical TKE models. 
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Table 3-1 

TYPICAL TKE MODELS FOR FREE SHEAR LAYERS 


1 . 


❖ 

The Prandtl Energy Model 


(Ref. 3-12) 


- <VF 1 VF2> = 0.08 



[o.625 XG2 


9VM1 

9X2 


(3.65) 


where XG2 is a lateral distance between the dimensionless 0.1 and 0.9 
velocity levels in the free shear layer. 
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(3.66) 


where 


e = 0.08 p (K/2) 2 /E (3.67) 

❖ 

2. Energy Dissipation Model (Ref. 3-12) 
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DE 
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f e 
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IS) + 2 - 86 l e 

&m 2 - ^(¥) 2 

(3.70) 


3. Dynamic Equations with algebraic length scales, Ref. 3-9. 


4. Dynamic Equations with another partial differential equation for the decay 
of turbulent kinetic energy (Ref. 3-12) 


Axisymmetric versions of 1 and 2 are also reported. 


3-29 

LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC-HREC D306102 


Section 4 

EDDY VISCOSITY MODELS 
4.1 FREE SHEAR LAYERS 

The prediction of shear layer properties for reacting flows requires 
empirical eddy viscosity coefficient models be used in order to remain within 
practical limits. Even then only parallel flows can be described. In the 
process of solving problems such as plume afterburning and jet engine or 
rocket plume dissipation, a working principle has evolved for selecting among 
the many eddy viscosity models the ones to be applied to particular types of 
shear layers. 

An eddy viscosity model, which would treat all the mixing problems 
encountered, is desirable. This goal has not been achieved. If a single 
model must be used for all problems, a choice could of course be made, 
but in certain regions or for certain types of flows, the knowledge that a 
different eddy viscosity model performs more realistically will invariably 
mean that the more realistic model will be used to treat that particular re- 
gion or problem. For that reason discussion of empirical eddy viscosity 
models should be based upon the kinds of problems a model will be used to 
solve . 


The eddy viscosity models described will be of the following form 


P € 



(4.1) 


where p is the density, € is the turbulent kinematic eddy viscosity, C is an 
empirically determined constant, L is a length characteristic of the mixing 
and the fjpuj implies that the eddy viscosity models are functions of the mass 
flux. Laminar mixing can be represented by a special form of this mixing 
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model; i.e., 

pe = Cj (4.2) 

where C at most can be a function of temperature and composition. 

The choice of which form, turbulent or laminar, is determined by the 
particular state of the flow being calculated. 


4.1.1 Planar Mixing 


An eddy viscosity model of the Prandtl type is used to describe the 
shear layer which develops between two planar streams, 1 and 2. 


P« 


0.0011 At 


m 


Pl U l 


-P 2 U 2 


(4.3) 


where At^ is the thickness of the mixing region. The outer edge of the mixing 
region is defined as the position where the mass flux ratio of successive 
stream function positions, (p.U. )/p. , .U. , . ) is 0.999. The inner edge is defined 

IX 1 I x 1 r 1 ^ 2 

similarly using the ratio (p i+1 lL + j )/p^ U-. The units on p are slug/ft ; on ( , ft / 
secjonAt^, ft; and on U, ft/sec^. These units result in a pe product with units 
of slug/ft/sec. These can be converted to the more common units of lb^-sec/ft 
by dividing the entire equation by the gravitational constant, gc where 


g c 


1 slug -ft/sec' 
1 lb. 


(4.4) 


Note that on occasions eddy viscosity models are listed which have similar 
form but greatly different constant values. The most probable explanation 
is that a different set of units is used rather than a different magnitude of 
eddy viscosity. 


Many values of the constant in Eq. (4.3) have appeared over the years. 
These relatively small variations have been attributed to better matches of 
particular sets of experiments similar to those of interest to the user. The 
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constant value 0.0011 in Eq, (4.3) was obtained from studies of mixing and 
combusting two-dimensional flows (Ref. 4-1). These are the types of flows 
of greatest interest to this particular user, hence the 0.0011 value is rather 
arbitrarily listed. 

In studies on chemical combustion lasers involving two-dimensional 
mixing and combustion, a parametric variation on eddy viscosity models was 
run by changing the constant values up to an order of magnitude (Ref. 4-2). 

The end result of the mixing and combustion (in this case energy in the in- 
verted population state) varied by only 15% of the peak potential power of the 
laser. The amount of change for any specific problem is obviously dependent 
on the features of the problem. In studying a wide variety of problems, how- 
ever, experience has shown that problem modeling has been more important 
than minor variations in the value of the constant in the eddy viscosity repre- 
sentation. The skill, in setting up the problem (within the often rather severe 
restrictions of the mixing analyses) to best describe the real flow process 
which occurs, can have a significant effect on the results of the mixing study. 
Concisely, one should recognize that a complicated process is being approxi- 
mated simply, and should be prepared to accept the results as good indications 
of those processes for which they were developed. It should nevertheless be 
recognized that a certain degree of the user's judgment is required to interpret 
those results properly. 

4.1.2 Axisymmetric Mixing 

The shear layer which develops as an axisymmetric jet mixes with its 
surroundings, whether quiescent or co-flowing, is calculated by using the 
combination of two-eddy viscosity models. The first model is applied in the 
initial region of the jet where the shear layer is developing. The second 
eddy viscosity model is used in the developed region of the flow where no 
potential jet core remains. This approach reduces the transition zone to 
zero and directly connects the two regions. 
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The eddy viscosity model recommended for the initial region has the 
same form as that suggested for the two-dimensional shear layer but has a 
different value for the constant (Ref. 4-3). 


PC 


= 0.01 At 


m 


P U 

I oo oo 


P «L U <L 


(4.5) 


This eddy viscosity model is used until the jet centerline properties 
have been changed by the mixing (i.e., the potential core region is gone). 

The entire mass flux of the jet is then involved in the shear layer. An eddy 
viscosity model more suited to thick, developed shear layers is applied for 
the remainder of the shear layer calculations. A Ferri eddy viscosity model 
is used to predict this regime. The Ferri model has the form 


pc 


0.018 


r l/2 



- p «, u «, 


(4.6) 


where the T ^/z va ^ ue t ^ le ra dial location in feet at which the local pU product 
is 


pU = l/2(p <c U o+ p 1J .U (S .) 


(4.7) 


Recalling Eq. (4.1) and Eq. (4.5), the form of Eq. (4.6) can be altered 
to make the changeover of eddy viscosity models appear less abrupt. Re- 
writing Eq. (4.6) in the form of Eq. (4.1), 


p€ = 0.01 (1.8 r l/2 ) jp^U^ - Pff U 


00 00 




(4.8) 


A balanced axisymmetric jet mixing region is sketched in Fig. 4-1 . 

At the end of the potential core region the At m for a balanced jet (to which the 
equations describing the mixing, restrict in the strict sense all mixing calcu- 
lations) is approximately two radii thick. From Eq. (4.8) the (1.8 r^^) ^ actor 
makes the switchover of eddy viscosity models appear to be more consistent. 
The eddy viscosity model used in the initial region is most applicable for 
relatively thin shear layers. Questions arise as one continues to use this 
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Fig. 4-1 - Balanced Axisymmetric Jet Mixing Region 

model until either the entire potential core region has disappeared or the model 
predicts a viscosity level equal to the Ferri model for the developed region. 
Values in the vicinity of the transition from the initial region to the developed 
region should be judiciously considered because of the assumption of zero 
transition length. If the area of interest lies in the transition zone, some 
error estimates should be placed on the answers. 

Any shear layer predictions close to the starting of the calculation should 
be viewed with caution. The discrimination of the edges of the mixing region 
are greatly influenced by the grid size in the first few steps of the calculation. 

If a simple case of a uniform jet and external flow is being considered, the 
thickness of the mixing zone is forced to be (at least) one grid size thick. This 
mixing zone thickness is unfortunately a function of the numerical scheme of 
solution and not of the physical processes. As the calculation proceeds down- 
stream, the viscous action replaces the grid size as the driving function of the 
thickness of the mixing zone. 

As a working rule, the mixing calculation results are of the greatest 
value in that ambiguous region, sufficiently far downstream of the initial 
region and not so far downstream to have lost all influence of the input 
conditions. 
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Since the recommended models were developed by matching experi- 
mental data, the models have been used with success to treat even com- 
busting flows. Because of the difficulty of realistically relating turbulent 
shear stress energy levels and the combustion processes, the empirical 
eddy viscosity model will probably continue to be used for a considerable 
period when combusting mixing processes are calculated. 

4.2 BOUNDARY LAYERS 


Boundary layer flows have been calculated using eddy viscosity coeffi- 
cients. One of the more successful models, Cebeci's extended eddy viscosity 
model (Ref. 4-4), uses a two-layer representation of the eddy viscosity. In 
the inner region, closer to the wall, the eddy viscosity is based on Prandtl's 
mixing length theory, as modified by Van Driest, to account for the damping 
effect of the wall, and as extended by Cebeci to include wall mass transfer, 
compressibility and pressure gradient effects. The eddy viscosity in the 
inner region is given by 


€ . = r 

i 


au 


dY 


(4.9) 


where the mixing length, l, is 


f = 0.4 Y [l - exp(-Y/A)] 


(4.10) 


the Van Driest damping factor A is defined as 

A = 26_u 

(T w P) l/Z N 


(4.11) 


and the factor N which accounts for pressure gradient and mass transfer 
effects is given by 


N 


2 


dP jl 1 
d X(pV) w T w 


- exp 


11.8 (pV) U 11.8 (pV) p 

V* ' w V- ' w “ w 

+ exp 


^w T w )1/Z f 


(p T ) 1//2 P 
p w w ^ 


(4.12) 
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If (pV) w equals zero N is calculated as 


N 2 = 1 + 11.8 — (p f 1 
dX 


/ 

\3 


Y 

^w ( 

VST 



w 


1/2 


(4.13) 


In the outer, wake-like, portion of the boundary layer, Clauser's form of the 
eddy viscosity, modified to include an intermittency factor, is used. The 


outer eddy viscosity is given by 
€ = 0.0168 U 


/(•■it) 


dY 


1 + 5.5 (Y/6) 


1 " 


LO 


(4.14) 


where the term in brackets is an approximation to Klebanoff's error function 
intermittency relationship. 

The eddy viscosity for the inner region is used from the wall outward 
until the height at which £ q = ^ is reached. From that point to the boundary 
layer edge the outer expression for eddy viscosity is utilized. 
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Section 5 

CONCLUSIONS AND RECOMMENDATIONS 


The work reviewed in this report leads one to the following conclusions 


1. Use of the dynamic equations for single-point correlations 
and of a transport equation for the mixing length as defined 
by a two-point correlation is a major improvement in tur- 
bulence models. 

2. Insufficient results are presently available to exploit fully 
the models mentioned in Conclusion (1). 

3. Present simple turbulent kinetic-energy models offer no 
significant advantage over eddy -viscosity models for free- 
shear layers. 

4. The only possible advantages which TKE models could 
offer for free shear layers is to predict a realistic lateral 
and streamwise variation of mixing rates. Such variations 
could be used empirically with eddy -viscosity models, if 
good experimental shear layer descriptions existed. These 
statements are particularly true for combusting jet flows. 

5. As a first step, intermittence factors which have been 
evaluated for boundary layer flow, could be used directly 
in eddy viscosity models. 


The following recommendations are offered. 


1. Those interested in turbulent mixing calculations should 
stay abreast of TKE model development. 

2. Eddy viscosity models which vary laterally should be de- 
veloped and their behavior investigated, particularly for 
combusting flows. 
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